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ABSTRACT 

We investigate the accuracy of the parametric recovery of the line-of-sight velocity distribution (LOS VD) of the 
stars in a galaxy, while working in pixel space. Problems appear when the data have a low signal-to-noise ratio, 
or the observed LOSVD is not well sampled by the data. We propose a simple solution based on maximum 
penalized likelihood and we apply it to the common situation in which the LOSVD is described by a Gauss- 
Hermite series. We compare different techniques by extracting the stellar kinematics from observations of the 
barred lenticular galaxy NGC 3384 obtained with the SAURON integral-field spectrograph. 
Subject headings: galaxies: individual (NGC 3384) - galaxies: kinematics and dynamics - line: profiles 
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1. INTRODUCTION 

The dynamics of stars in a galaxy is uniquely defined by 
the distribution function and the potential in which the stars 
move. From the many stable equilibrium configurations for 
collisionless systems which can in principle be constructed, 
only some of them are actually observed in our Universe. This 
constitutes a 'fossil record' of galaxy formation. It is still 
not clear what observable constraints are needed to recover 
both the distribution function and the gravitational potential 
of a stellar system (but see, e.g., Dejonghe & Merritt 1992, 
for spherical geometry), but simple dimensionality arguments 
imply that this should not be possible without the knowledge 
of the full line-of-sight velocity distributions (LOSVD) at all 
spatial positions on the galaxy image on the sky. It is therefore 
useful to explore how the LOSVD can be recovered from the 
observations. 

Considering galaxies as pure stellar systems, the spectrum 
observed at a certain sky position is a (luminosity-weighted) 
sum of individual stellar spectra redshifted according to their 
line-of-sight velocities. If one makes the assumption that the 
spectrum of all stars is given by a single template, then it 
simply reduces to the convolution between that spectrum and 
the LOSVD, which can then be retrieved by solving the in- 
verse problem, i.e., deconvolving the spectra using the tem- 
plate. Deconvolution is an intrinsically ill-conditioned prob- 
lem, which amplifies noise and measurement errors. Due to 
the special care required in the inversion many techniques 
have been developed in the last 30 years to recover the 
LOSVD from the data. The evolution in the methods has 
been mainly driven by the gradual improvement in the obser- 
vational techniques and the signal-to-noise ratio (S/N) of the 
data, and by the steady increase in the available computational 
speed. 

Early methods were mostly using Fourier based techniques, 
which allowed the LOSVD to be recovered quickly, from a 
deconvolution process, and in some cases included techniques 
to reduce the effect of template mismatch (Simkin 1974; Sar- 
gent et al. 1977; Tonry & Davis 1979; Franx & Illingworth 
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1988; Bender 1990; Statler 1995). More recently methods 
have shifted towards fitting the LOSVD directly in pixel space 
(Rix & White 1992; Kuijken & Merrifield 1993; van der 
Marel 1994; Saha & Williams 1994; Merritt 1997; Gebhardt 
et al. 2000; Kelson et al. 2000). The reasons for this are that (i) 
in pixel space it becomes easy to exclude gas emission lines 
or bad pixels from the fit, and take continuum matching di- 
rectly into account; (ii) current computers can accommodate 
the larger computational cost involved; (iii) the availability of 
libraries with high spectral resolution stellar and galaxy spec- 
tra allows the template to be carefully matched to the observed 
galaxy spectrum (e.g., Emsellem et al. 2004). See de Bruyne 
et al. (2003) for a more detailed historical overview of the 
various methods. 

The different techniques can be further subdivided accord- 
ing to whether the LOSVD is derived in a non-parametric way 
(in practice computed on a small set of discrete values) or 
parametrically, as a function of a limited number of param- 
eters. In the latter case the Gauss-Hermite parametrization 
by van der Marel & Franx (1993) and Gerhard (1993) is es- 
sentially always adopted (but see, e.g., Zhao & Prada 1996, 
for an alternative). However even when the LOSVD is deter- 
mined non-parametrically, the Gauss-Hermite parameters are 
still generally used to present the result in an easily under- 
standable way. 

In this paper we study again the problem of recover- 
ing, while working in pixel space, a LOSVD described by 
the Gauss-Hermite parametrization. Compared to the non- 
parametric case, the process and the estimation of measure- 
ment errors are simplified. In Sec. 2 we describe the general 
problem. In Sec. 3 we discuss different approaches to the ex- 
traction, we find that special care has to be taken when the 
LOSVD is undersampled by the data or the S/N is low, and 
we present a solution based on the maximum penalized like- 
lihood formalism. In Sec. 4 we draw some conclusions. 

2. FORMULATION OF THE PROBLEM 

The parametric recovery of the LOSVD in pixel space starts 
with creating a model galaxy spectrum G mo d(x), by convolv- 
ing a template spectrum T(x) by a parametrized LOSVD. 
Both the object and the template spectra are rebinned in wave- 
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length to a linear scale in x = In A, while usually preserving 
the number of spectral pixels. The best-fitting parameters 
of the LOSVD are determined by minimizing the % 2 , which 
measures the agreement between the model and the observed 
galaxy spectrum G(x), over the set of N good pixels: 



x 2 =E^ 2 



where the residuals are defined as 

G"mod(Xtt) - 



G(x n ) 



(1) 



(2) 



AG(x n ) ' 

with AG(x„) the measurement error on G(x n ). 

More specifically the following model is adopted for the 
galaxy spectrum: 

K L 

Gmodto = ^tv*[5 * T k ](x)+J2biVi(x) w k > 0, (3) 



k=l 



1=0 



where 71 is in general a library of K galaxy or stellar tem- 
plates, B(x) = C(cx) is the broadening function, with £(v) the 
LOSVD, c is the speed of light and * denotes convolution. 
The Vi(x) are here chosen to be the Legendre polynomials of 
order I and account for low frequency differences in shape be- 
tween the galaxy and the templates. For each given C(v), the 
optimization of x 2 is a bounded-variables linear least-squares 
problem for the weights (w\,. . . ,WK,bo, ■ ■ ■ ,b£) which can be 
solved, e.g., with the specific BVLS algorithm by Lawson & 
Hanson (1995), or as a quadratic programming problem. Here 
we are interested in the determination of the parameters defin- 
ing C(v) and in what follows we will assume that the weights 
of Eq. (3) are always optimized in this way. Multiplicative 
polynomials can also be included in the fit (see Kelson et al. 
2000), without affecting the discussion that follows. 

Following van der Marel & Franx (1993) and Gerhard 
(1993) it has become standard to expand the LOSVD as a 
Gauss-Hermite series 



Av) = 



,-(1/2)/ 



l+^2h m H m (y) 



m=3 



(4) 



where y = (v-V)/a and the H m are the Hermite polynomials. 
With these definitions the minimization of the \ 2 m Eq. (1) 
is a nonlinear least-squares optimization problem for the M 
parameters (V, a, hi,... ,Hm)- Least-squares problems can be 
solved much more efficiently than general ones, by using spe- 
cific algorithms which require the user to provide the residuals 
r n of Eq. (2), to compute explicitly the Hessian matrix of the 
X 2 merit function (see, e.g., Press et al. 1992, §15.5). Here we 
will use the MINPACK 1 implementation (More et al. 1980) of 
the Levenberg-Marquardt method for nonlinear least-squares 
problems. 

3. DISCUSSION 

In this section we compare three different approaches to the 
determination of the best-fitting parameters of the LOSVD in 
Eq. (4). We explain the limitations of the different methods 
and finally select the last one as the optimal choice. We do not 
address here the template mismatch issue, which we assume 
is minimized by the choice of an optimal library of templates 
in Eq. (3) as in Emsellem et al. (2004). 

1 We used an IDL porting of the code by Craig B. Markwardt, available 
from the Web page http://cow.physics.wisc.edu/~craigm/idl/ 
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FIG. 1. — Double Gaussian representation (Eq. 5) of the LOSVD ob- 
served by Scorza & Bender (1995) at 17" along the major axis of the SO 
galaxy NGC 3115 (solid line). The two individual Gaussian components are 
shown with the gray lines, while the best fitting fourth order Gauss-Hermite 
parametrization is plotted with the dashed line. 



To explain the characteristics of the three methods, we will 
use each of them to extract a realistic but known LOSVD, 
observed by Scorza & Bender (1995) at r = 17" along the ma- 
jor axis of the SO galaxy NGC 3115 (their Fig. A.2). This 
LOSVD (Fig. 1) is representative of the one observed in a 
number of galaxies, and can be very well described by a dou- 
ble Gaussian parametrization 



F(v) = ^/,-exp 



7=1 



(y-Vjf 

2 a) 



(5) 



with parameters h = 0.041, 1 2 = 0.032, Vi = 48.7 km s~\ V 2 = 
-77.3 km s" 1 , a x = 70.0 km s" 1 and ct 2 = 130.0 km s" 1 (the 
LOSVD was shifted to zero mean velocity). 

The best fit to F(y) using a fourth order Gauss-Hermite 
series (Eq. 4) is obtained with parameters V = 0.0 kms -1 , 
cr = 1 14.8 km s~\ h 3 = -0.150 and h 4 = 0.036, while the best 
fitting Gaussian has V = -2.3 km s -1 and a = 1 18.9 km s" 1 . 

3.1. Fitting (V, a) first 

In an ideal situation where the LOSVD is perfectly sam- 
pled by the data, the (h 3 , . . . ,h M ) parameters of the LOSVD 
(Eq. 4) are essentially uncorrected to (V,a). Therefore one 
expects the best-fitting parameters to change little, irrespec- 
tive of whether (V, a) are fitted first or together with the other 
parameters. To lowest order the difference between the two 
approaches is given by Eq. (10) of van der Marel & Franx 
(1993). 

An advantage of fitting (V, a) first is that the Gauss-Hermite 
parameters that one obtains, coincide with the true Gauss- 
Hermite moments integrated over the LOSVD (see van der 
Marel & Franx 1993, for a detailed discussion). This also 
means that the value of the parameters (hi,,..., h M ) do not de- 
pend on the adopted number M of terms in the expansion. 
This precise relation between the LOSVD and the extracted 
kinematical parameters can be very useful for the dynamical 
modeling (e.g. Rix et al. 1997). 

Moreover when only a single template is adopted (K = 1 in 
Eq. 3), for each pair of (V, cr) values the best fitting value of 
all the remaining parameters (w\ , hj, , . . . , Hm, bo, . . . , b£) can be 
determined linearly. These ideas led van der Marel (1994) to 
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FIG. 2. — Logarithmically rebinned galaxy model spectrum (thin noisy line) 
convolved with the LOSVD of Fig. 1 and with added Poissonian noise at a 
S/N = 60. The dots represent the residual difference between this model and 
the best fitting template (thick gray line). The spectrum covers the wavelength 
range 4800-5380 A, with a pixel scale of 60 km s -1 , and includes prominent 
absorption features of H/3, Mgb and Fe5270. 



design an efficient method to fit the parameters of the LOSVD 
in pixel space, where (V, a) are fitted first and the remaining 
parameters are linearly expanded at the best fitting (V, a) so- 
lution. 

For extracting the stellar kinematics of the 72 galaxies of 
the SAURON survey (de Zeeuw et al. 2002) we needed to op- 
erate in pixel space to be able to deal with the contamination 
due to emission-line gas present in most of the observed ob- 
jects (Emsellem et al. 2004). The speed of van der Marel 
(1994) pixel fitting algorithm was an attractive feature, given 
the large number (<~ 200.000) of independent spectra from 
the survey and the need to compute accurate errors by Monte 
Carlo simulations. However the observed LOSVDs are not 
always well sampled by the SAURON pixel scale of dv « 60 
km s" 1 . The Nyquist critical frequency for h\ is <~ er/2. This 
means that undersampling becomes a problem for the deriva- 
tion of the first four Gauss-Hermite moments of the LOSVD 
when the observed velocity dispersion is of the order of ~ 2 
pixels (for SAURON 120 km s" 1 ). 

To test the accuracy of the recovery of the Gauss-Hermite 
moments using this method, when undersampling becomes 
significant, we first created a realistic model spectrum by fit- 
ting a (logarithmically rebinned) library of Vazdekis (1999) 
galaxy model spectra to the average SAURON spectrum of the 
barred lenticular galaxy NGC 3384. This spectrum was over- 
sampled by smoothly interpolating it on a 30 times finer spec- 
tral grid and was subsequently convolved with the LOSVD of 
Fig. 1 . It was integrated over the SAURON pixels, in the wave- 
length range 4800-5380 A, and Poissonian noise was finally 
added to represent the simulated galaxy spectrum G(x). The 
template T(x) was obtained by integrating the original over- 
sampled spectrum over the SAURON pixels. In Fig. 2 we show 
an example of such a spectrum, with noise added at S/N = 60. 

We subsequently tried to recover the true Gauss-Hermite 
moments, by fitting (V, a) first and expanding (fi3,...,fiM) lin- 
early at the best fitting (V, a) location. To check our results we 
computed the true moments ([13,114), by analytically integrat- 
ing over the LOSVD (Eq. 8 of van der Marel & Franx 1993), 
using (V,a) of the best fitting Gaussian computed in Sec. 3. 
We found hj, = -0.144 and h.4 = 0.029, which as expected are 




FIG. 3. — Fitting (V,a) first. Recovery of the stellar kinematics from a 
model galaxy spectrum convolved with the LOSVD of Fig. 1. The shape of 
the LOSVD was held fixed, while the velocity scale was varied so that the a m 
of the best fitting Gauss-Hermite series varied in the range 48-360 km s _1 , 
corresponding to 0.8-6 pixels. In these measurements V and a were fitted 
first (nonlinearly) and only then and I14 were fitted (linearly) with V and a 
fixed to the optimal values. In each panel the thick dotted line represents the 
true Gauss-Hermite moments, that the method is trying to recover, the thick 
solid line indicates the value of the best fitting Gauss-Hermite parameters, 
while the crosses represent the actually measured values, as a function of 
the input <j m , for 1000 different Monte Carlo realizations with S/N = 60. 
The thick noisy line shows the measured values for 100 realizations with a 
very high S/N = 600. The recovered moments start becoming biased when 
(Tin i; 240 (4 pixels), and this method becomes essentially insensitive to any 
deviation from a Gaussian when <Tj n < 120 km s _1 (2 pixels), at any S/N. 



very close to the best fitting parameters computed before. 

The result of this test, both for a S/N = 60 and S/N = 600, 
is shown in Fig. 3. Here the shape of the input LOSVD was 
held fixed, but the velocity scale was varied so that the <j\ n 
of the best fitting Gauss-Hermite series varied in the range 
48-360 km s~\ corresponding to 0.8-6 SAURON pixels, for 
1000 different Monte Carlo realizations. By construction, for 
(Tin = 1 14.8 km s" 1 the LOSVD reduces to the one presented 
in Fig. 1. The recovery of the true moments is strongly bi- 
ased when (Ti n ^5 240 km s" 1 (4 pixels). This can be under- 
stood by the fact that, when the LOSVD is not very well sam- 
pled, an asymmetry in the profile can be compensated, during 
the Gaussian fit, by a small V shift, while a symmetric de- 
viation from a Gaussian is nearly equivalent to a change of 
a. The bias in the recovery of the input parameters remains 
unchanged even for a noise-free model spectrum, although 
of course the scatter in the values disappears. And no im- 
provement is observed if the template is oversampled before 
convolving it with a well sampled LOSVD in Eq. (3). These 
results are representative of what we found with different re- 
alistic input LOSVDs. 

What was seen in the simulation is also found on real 
SAURON data. The approach of fitting (V, a) first, becomes in- 
sensitive to any deviation from a Gaussian LOSVD, when the 
observed velocity dispersion is low (< 120 km s" 1 , or 2 pix- 
els). As a practical example this method gives the misleading 
impression that the LOSVD in the nucleus of a low-dispersion 
object like NGC 3384, containing a fast-rotating disk-like 
structure, is consistent with being essentially symmetric, as, 
e.g., the central LOSVD of the high-dispersion non-rotating 
giant elliptical galaxy M 87. This apparent similarity of 
the overall shape of LOSVD in the two different galaxies is 
an artefact of the extraction method, and we will show that 
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FIG. 4. — Same as in Fig. 3, but in this plot the four parameters (V, a, (13,114) 
were all fitted simultaneously (nonlinearly). Now the measurements are sup- 
posed to reproduce the best-fitting Gauss-Hermite parameters, which are in- 
dicated by the thick solid line. The measurements are generally unbiased, but 
no reliable and h$ measurements can be obtained at S /N = 60 for <x in < 120 
km s - ' (2 pixels). At low a- m the scatter in V and a increases due to their 
correlation with I13 and A4 respectively. Moreover the measured a tends to 
be systematically lower than the true value, while A4 is correspondingly too 
high. In the case of very high S/N the strong asymmetry of the LOS VD can 
be recovered, with a modest bias, down to the smallest u m values. 



it is still possible to recover, with a different method, the 
strong asymmetry of the LOSVD from the SAURON data of 
NGC 3384, even at low dispersion. 

3.2. Fitting all parameters simultaneously 

The bias present in the recovery of the Gauss-Hermite mo- 
ments of an undersampled LOSVD, in the previous section, 
was mainly due to the fact that (V, a) were not fitted simulta- 
neously with the other LOSVD parameters. The fit should 
then become generally unbiased when all (V,<7,/j 3 , . . . ,h M ) 
parameters are varied simultaneously to optimize the fit. In 
this way one should recover the best fitting parameters of the 
Gauss-Hermite series, which do not precisely coincide with 
the true Gauss-Hermite moments of the LOSVD. By defini- 
tion this method will provide the lowest x 2 m the fit, and for 
this reason it was also originally chosen by van der Marel & 
Franx(1993). 

In Fig. 4 we repeated the experiment of the previous sec- 
tion while fitting all the above nonlinear parameters together. 
As expected the solution is now almost unbiased and the mea- 
surements tend to be spread around the true input values, for 
a wider a m range. However, when <T ln < 120 km s" 1 (2 pix- 
els) and the LOSVD is correspondingly not well sampled, the 
spectrum does not contain enough information to constrain 
all the free parameters and the scatter in the recovered values 
increases dramatically. 

Moreover, at low a m , the a and h\ values are not unbi- 
ased any more, not being symmetrically distributed around 
the known input values. The reason for this asymmetry can 
be understood by looking at the shape of the x 2 contours, 
which measure the agreement between the input spectrum 
G(x) and its best fitting model G mo d(x). In Fig. 5 we plot the 
Ax 2 = X 2 ~ Xmin contours for fits obtained by convolving G(x) 
with the LOSVD of Fig. 1, while keeping V and hj, fixed to 
the known best fitting values (Sec. 3) and varying a and hi, 
of G mo d(x). The narrow curved valley of nearly constant x 2 , 
for decreasing a and increasing hi from the location of the 
minimum Xmin< 15 due to the strong correlation between these 
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FIG. 5. — Contours of the A\ 2 = X 2 ~Xmin which measure the agreement 
between a simulated noiseless input galaxy spectrum G{x) and its best fitting 
model GmodOt). As in Fig. 4 G(x) was convolved with the LOSVD of Fig. 1. 
The Ax 2 is plotted as a function of a and A4 of the LOSVD of G mo i(x), 
while V and A 3 are kept fixed to the known best fitting values (Sec. 3). The 
plus sign indicates the location of the known best fitting (o-,/i 4 ). The three 
lowest Ax 2 levels correspond to confidence levels of 1, 2 and 3cr (thick line) 
respectively, for a S/N = 60, while other levels are separated by factors of 2. 
Note the narrow valley of nearly constant x 2 , for decreasing a and increasing 
/i 4 , which is due to the undersampling of the LOSVD. At this a m a very high 
S/N is needed for the fit to converge to the true minimum. 



two parameters, which is caused by the undersampling, and 
produces the effects seen in Fig. 4. 

When the errors in the measured parameters are large, one 
has to make a scatter versus bias tradeoff decision. Usually 
when the data are unable to tightly constrain a large number 
of parameters one prefers to retain in the fit only the parame- 
ters that are required to significantly decrease the x 2 - For this 
reason, while fitting a parametric LOSVD, it is common prac- 
tice to reduce the fit to (V, a) at low S /N: although the fit may 
be biased by the lack of freedom in the model, an additional 
flexibility would not lead to meaningful measurements, and 
structures and trends in the (V, a) values can still be detected 
due to the decreased scatter. 

Another problem of the large scatter that is present when fit- 
ting all parameters together is specific to the two-dimensional 
(2D) kinematical measurements obtained with an integral- 
field spectrograph. When presenting kinematics maps there 
is no easy way to also attach errors to the displayed values, 
and it becomes difficult to estimate from the maps when some 
structure is real and when it is due to noise. For this reason 
it would be preferable to display on the kinematics maps only 
the features that are statistically significant. 

3.3. Penalized pixel fitting 

The LOSVD of galaxies is generally well reproduced by 
a Gaussian (e.g., Bender, Saglia, & Gerhard 1994). For this 
reason one would like to use a technique to fit the LOSVD in 
which the solution is free to reproduce the details of the actual 
profile when the S/N is high, but where the solution tends to a 
Gaussian shape in case the S /N is low. This was done for the 
case the LOSVD is described by a non-parametric function, 
by using the maximum penalized likelihood formalism (e.g., 
Merritt 1997), and we refer to that paper for details. Here we 
apply the formalism to the case where the LOSVD is para- 
metrically expanded as a Gauss-Hermite series. We will show 
that this leads to a much simpler and faster implementation 
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than for the general case. 

The idea is to fit the parameters (V, a, hi,... ,Hm) simulta- 
neously as in Sec. 3.2, but to add an adjustable penalty term 
to the x 2 , to bias the solution towards a Gaussian shape, when 
the higher moments are unconstrained by the data, so that the 
penalized x 2 becomes: 

■X 2 + aV. 



4- 



(6) 



A natural form for the penalty function V is given by the in- 
tegrated square deviation of the line profile £(v) from its best 
fitting Gaussian Q(v): 



V 2 = 



Q£(v)-g(v)] 2 dv 



(7) 



This penalty does not suppress noisy solutions, which are al- 
ready excluded by the use of a low order parametric expansion 
for C{v). It was shown by van der Marel & Franx (1993) that 
in the case where C(v) has the form of Eq. (4) then Eq. (7) is 
well approximated by: 
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FIG. 6. — Penalized pixel fitting. Same as in Fig. 4 but using the new 
penalized parametric pixel fitting method with A = 0.7. At S/N = 60 the bias 
on the parameters is now significant only when a m goes below about 120 
km s~' (2 pixels). But more importantly the fit tends to converge to the true 
solution even for low <r m at the higher S/N = 600. 



m=3 



(8) 



In principle one could then define the penalty asV = T> z and 
optimize X 2 , (Eq. 6). In practice however this is not desirable 
for two reasons: 

• as explained in Sec. 2, it is computationally much more 
efficient to minimize the residuals r n (Eq. 2) instead of 
explicitly compute the x 2 ; 

• one needs a way to automatically adjust the penalty fac- 
tor a according to the x 2 of the observed fit. 

We found that a simple and effective solution to these prob- 
lems consists of using the following perturbed residuals as 
input to the nonlinear least-squares optimizer: 

r' n = r n + Xa(r)V, (9) 

where the variance is defined as 



(10) 



The qualitative interpretation of this formula is that a devia- 
tion V of the LOSVD from a Gaussian shape will be accepted 
as an improvement of the fit, only if it is able to correspond- 
ingly decrease the scatter cr(r) by an amount related to V. To 
quantify one can compute the objective function of the fit 



N N 

X 2 = Y/l + 2\a(r)VY /n + N [Xa(r)V] 2 



n=\ 



n=l 



(ID 



The sum of the residuals in the second term is zero by con- 
struction, due to the fact that the weights are optimized for 
given C{v) (Sec. 2; this is required for this form of perturba- 
tion to work). Considering the definition of variance (Eq. 10) 
one can finally write 

X 2 (l + A 2 2? 2 ), 



Xp 



(12) 



which is of the desired form (6) with a = X x automati- 
cally scaled according to the x 2 of the fit. In practice a(r) 
in Eq. (9) is computed using a robust biweight estimator 
(Hoaglin, Mosteller, & Tukey 1983), so that Eq. (12) is only 



valid as an approximation. One can see from this formula 
that, with A = 1, a deviation V = 10% of the LOSVD from a 
Gaussian (e.g., A3 = 0.1 and h m = 0) requires a corresponding 
decrease in the scatter cr(r) of the unperturbed residuals by 
more than (1+V 2 ) 1 / 2 = 0.5% to be accepted by the optimiza- 
tion routine as an improvement of the fit. 

A statistical interpretation of Eq. (12) can be obtained by 
considering that for a good fit x 2 ~ N. For a given A a devi- 
ation T> from a Gaussian needs to decrease x 2 by more than 
Ax 2 ~ N\ 2 T> 2 to be accepted, and this variation can be as- 
sociated to a specific confidence level (see, e.g., Press et al. 
1 992, § 1 5 .6) at which a Gaussian shape is excluded. Although 
this formula may serve as a guideline, we believe it is safer to 
test a choice of A using simulations as the ones presented. 

A test of the application of Eq. (9) to the same model spec- 
tra of Sec. 3.2 is shown in Fig. 6. For the same S/N = 60 as 
before, using A = 0.7, the measurements are now essentially 
unbiased when a m > 120 km s" 1 (2 pixels). For lower input 
dispersion the LOSVD converges towards a Gaussian and, as 
expected, the bias becomes similar to Fig. 3. However the 
crucial difference with the previous case consists of the fact 
that for higher S/N = 600 the solution tends to converge to 
the known best fitting values, due to the fact that the penalty 
scales with the x 2 - 

We mention here that alternative forms to Eq. (9) are pos- 
sible, with the requirement that the objective function has the 
form (6). One possibility is to include a multiplicative pertur- 
bation of the residuals: 



r' n = r n (l + XV 2 ). 
In this case the objective function becomes 



Xp 



X 2 (1+2AP 2 + A 2 P 4 ), 



(13) 



(14) 



and if one neglects the small term containing T> 4 , this becomes 
again of the form (6), with a = 2 Ax 2 - We verified that in prac- 
tice this alternative form of perturbing the residuals produces 
virtually the same results as using Eq. (9), if the same a pa- 
rameter is adopted. We prefer the perturbation (9) given its 
robustness against outliers, due to the use of the biweight in 
the computation of <r(r). The form (13) may be useful in im- 
plementations where the average residuals are not zero for a 
given C(v). 
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FIG. 7. — Comparison between the Voronoi 2D-binned (Cappellari & Copin 2003) SAURON stellar kinematics of the barred lenticular galaxy NGC 3384 
extracted with the different methods explored in this paper. Top panels: kinematics extracted by first fitting (V,cr) (nonlinearly) and then expanding the (/13 . A4) 
parameters at the optimal (V, cr) location. The Gauss-Hermite parameters are everywhere significantly suppressed by this method. Middle panels: kinematics 
obtained by fitting (V,a,h^,h^) simultaneously, without penalty. Due to the (T-/14 correlation (Fig. 5), caused by the undersampling of the LOSVD, the a is 
noisier and depressed, while /14 tends to be strongly positive. The /13 values fluctuates in the outer parts between large positive and negative values. Bottom 
panels: the kinematics measured with the new penalized PXF method (with A = 0.7). This overcomes the problems of the two previous approaches. Only the 
statistically significant non-Gaussian features are preserved, otherwise the solution is smoothly reduced to a Gaussian. Contours of the reconstructed surface 
brightness are superimposed (in 1 mag arcsec -2 steps). 



In Fig. 7 we compare the three different approaches to the 
recovery of the LOSVD by using them to extract the stel- 
lar kinematics from SAURON (Bacon et al. 2001) integral- 
field spectroscopic observations of the galaxy NGC 3384 (see 
de Zeeuw et al. 2002, for a description of the observations), 
which has a rather low velocity dispersion over the whole ob- 
served field. The data cube was spatially binned to a minimum 
S/N = 60 using the Voronoi 2D-binning method by Cappellari 
& Copin (2003). The same differences that we observed using 
synthetic model spectra can also be seen from real data: 

• The Gauss-Hermite parameters are everywhere signif- 
icantly suppressed when fitting V and a first (top pan- 
els); 

• the simultaneous fit of all parameters (middle panels) 
tends to produce a velocity dispersion map which is 
noisier and has smaller values, while tends to be 
strongly positive. The I13 values fluctuate in the outer 
parts between large positive and negative values; 

• the use of a penalty function (bottom panels) overcomes 
the problems of the two previous approaches. Only the 
statistically significant non-Gaussian features are pre- 
served, otherwise the solution is smoothly reduced to a 
Gaussian. 

We also applied the three methods to the extraction of the 
kinematics of all the 48 early-type galaxies of the SAURON 



survey (de Zeeuw et al. 2002). As expected from the simula- 
tions, the three techniques gave very different results only for 
galaxies like NGC 3384, which have a low velocity disper- 
sion. The kinematics extracted from SAURON observations 
of galaxies with dispersion > 180 km s" 1 over the whole ob- 
served field provided very similar result with the three meth- 
ods (though by construction not identical). 

3.4. Errors on the fitted parameters 

The availability of fast computers makes Monte Carlo sim- 
ulations the preferred way to estimate measurements errors 
for LOSVD extraction methods. This consists of repeating 
the full measurement process for a large number of different 
realization of the data, obtained by adding noise to the origi- 
nal spectra. 

The method discussed in this paper, like many other avail- 
able methods for extracting the LOSVD, makes use of some 
form of filtering or penalization of the solution, thus biasing 
the results to suppress the noise. It may be worth empha- 
sizing that, when computing Monte Carlo errors, one has to 
correct for the bias in the measurements, introduced by the 
noise suppression mechanism. Failure to correct for the bias 
in the kinematics can provide a significant underestimation of 
the confidence intervals of parameters determined by fitting 
dynamical models to the data. 

One way to obtain proper errors from penalized methods 
consists of correcting the error estimates by increasing the 
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percentile intervals by an amount given by the bias (e.g., 
Beers, Flynn, & Gebhardt 1990). In real measurements one 
does not know the actual bias (otherwise the measurements 
would have been corrected for it), so an estimate of the maxi- 
mum expected bias may be used instead. This requires one to 
perform Monte Carlo simulations of the kinematic extraction 
to determine the maximum bias at different S /N levels and 
ai n values. 

A simpler order of magnitude estimate of the errors on the 
parameters may be obtained by noting that Fig. 4 provides 
a more realistic representation of the scatter in the measure- 
ments than Fig. 6, at low a m values. In the case of the pe- 
nalized pixel fitting routine described in the previous section 
one can then obtain the errors by setting the parameter A in 
Eq. (9) to a very small value before running the Monte Carlo 
simulations. This simple approach will generally provide a 
conservative estimate of the actual errors. 

3.5. The algorithm 

To summarize the discussion of the previous sections, the 
suggested algorithm for recovering the LOS VD parametrized 
as a Gauss-Hermite series is the following: 

1. start with an initial guess for the parameters (V, a), 
while setting the initial Gauss-Hermite parameters 
h 3 ,...,h M = 0; 

2. solve the subproblem of Eq. (3) for the weights 

(w 1 ,---,w K ,b ,...,b L ); 

3. compute the residuals r n from the fit using Eq. (2); 

4. perturb the residuals as in Eq. (9) to get r' n ; 

5. feed the perturbed residuals r' n into a nonlinear least- 
squares optimization routine and iterate the procedure 
from step (2), to fit for the parameters (V, a, hj, , . . . , h M )- 
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3.6. Availability 

An IDL routine implementing the algorithm described 
in this work is made available from the Web address 
http://www.strw.leidenuniv.nl/^mcappell/idl/. 

4. CONCLUSIONS 

In this paper we addressed the problem of extracting the 
LOSVD of the stars, parametrized using a Gauss-Hermite 
series, from observed galaxy spectra. Using a method that 
works directly in pixel space, we compared different tech- 
niques by applying them to the same problem, with both sim- 
ulated and real data, and we showed that one has to pay spe- 
cial attention to the extraction, when the LOSVD is not well 
sampled by the data, or when the S /N is low. 

We proposed that in these situations one should apply the 
maximum penalized likelihood formalism to extract as much 
information as possible from the spectra, while suppressing 
the noise in the solution. We demonstrated that this leads to 
a fast and simple algorithm. We also discussed how Monte 
Carlo errors should be properly estimated. This penalized 
pixel fitting method is particularly useful to extract the kine- 
matics from integral-field spectroscopic data, as for 2D maps 
there is no standard way to visualize errors and one wants to 
be able to show only the non-Gaussian features that are statis- 
tically significant. A routine that implements the ideas of this 
paper is made publicly available. 
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